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ABSTRACT 

We put forward and test a simple description of multi-point propagators (MP), 
which serve as building-blocks to calculate the nonlinear matter power spectrum. 
On large scales these propagators reduce to the well-known kernels in standard 
perturbation theory, while at smaller scales they are suppresed due to nonlinear 
couplings. Through extensive testing with numerical simulations we find that this 
decay is characterized by the same damping scale for both two and three-point 
propagators. In turn this transition can be well modeled with resummation results 
that exponentiate one-loop computations. For the first time, we measure the four 
components of the non-linear (two-point) propagator using dedicated simulations 
started from two independent random Gaussian fields for positions and velocities, 
verifying in detail the fundamentals of propagator resummation. 

We use these results to develop an implementation of the MP-expansion for the 
nonlinear power spectrum that only requires seconds to evaluate at BAO scales. To 
test it we construct six suites of large numerical simulations with different cosmologies. 
From these and LasDamas runs we show that the nonlinear power spectrum can be 
described at the < 2% level at BAO scales for redshifts in the range [0 — 2.5]. We 
make a public release of the MPTbreeze code with the hope that it can be useful to 
the community. 

Key words: cosmological perturbation theory - cosmological parameters - baryon 
acoustic oscillations - large-scale structure of the universe 



1 INTRODUCTION 

Ongoing and future galaxy redshift surveys will render the 
large scale structure of the Universe with unprecedented de- 
tail thanks to a combination of redshift depth and large 
survey area. Among them are the Sloan Digital Sky Sur- 
vey IlQ the WiggleZ survejj^] the Dark Energy Surve^J^] 
the Physics of the Accelerating Universe collaboratiorj^] and 
ES A/Euclid survey^] The main driver underneath this ef- 
fort is to seed light into the present cosmic acceleration. Var- 
ious probes exist that connect different statistical aspects of 
galaxies properties to cosmological parameters, in particular 
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to those related to acceleration, such as the Baryon Acous- 
tic Oscillations (BAO), Redshift Space Distortions (RSD) 
or Weak Lensing (WL). However in order to maximize the 
scientific outcome from this data we need to put forward 
precise theoretical and/or numerical predictions, for exam- 
ple for two-point statistics. The difficulty arise because the 
most rewarding range of scales lie in the nonlinear regime of 
structure formation. 

For BAO (and RSD) these nonlinearities, due in large 
part to gravitational instabilities, are not strong (Eisen- 



stein ct al. 2007 , Croccc & Scoccimarro 2008 , Angulo ct al 



2008 1 . Hence the problem can be addressed using pertur- 



bative schemes of the equations of motion in addition to 
numerical N-body simulations. This is then the main goal of 
this paper, to propose an implementation of a (resummed) 
perturbative expansion for the matter power spectrum that 
is both accurate (percent level) and practical (few seconds 
of evaluation) on large BAO scales. And to have it tested as 
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Run Q m Of, h erg n s L box (h 1 Mpc) N rU ns 

FID 0.27 0.04 0.7 0.9 1 1280 50 

tilt 0.9 1250 4 

WMAP3 0.2383 0.0418 0.73 0.74 0.95 1250 4 

Low-Q m 0.20 1250 4 

Mid-o-8 0.8 1250 4 

Low-cr 8 0.7 1250 4 

LasDamas 0.25 0.04 0.7 0.8 1 2400 4 



Table 1. Our suite of N-body simulations spanning different cosmological models. The top entry corresponds to our largest ensemble 
used as a benchmark for testing different components of the model. The last entry refers in particular to four Oriana runs, the largest 
box-size within the LasDamas simulations. Null entries indicate the same value as the FID run. 



much as possible against numerical simulations of different 
cosmological models. 

The process of nonlinear structure formation can be 
very well traced by gravitational N-body codes once certain 
requirements are met (e.g. Heitmann et al. (20101 and refer- 
ences therein). These involve, among others, good mass res- 
olution (i.e. particle load), fine time stepping, high starting 
redshift, large box-size (to include long wavelength modes), 
etc. A percent level estimate of the power spectrum puts 
strong constrains in these parameters which generally slow 
down the numerical solver. In addition it is numerically 
very expensive to explore cosmological parameter space with 
large high resolution simulations, although efforts in this di- 
rection are currently ongoing ( Heitmann et al.|2 009). 

In turn, cosmological perturbation theory (PT) is a 
well defined formalism that can be applied without extra 
cost to any LCDM model (and even beyond) but leads 
to poorly convergent results (see Bernardeau et al. ( 2002 1 



for a review). One step beyond this issue was put forward 
by |Crocce fe Scoccimarro| ( |2006a|b] ) through a systematic 
re-organization of the perturbative series so called Renor- 
malized Perturbation Theory. This led to a better behaved 
expansion and more robust results ( jCrocce fe Scoccimarro| 
2008). A number of similar studies with alternative methods 



to resum the PT expansion quickly followed, e.g . |Matarrese 
fc Pietroni| (|2007|>; |Taruya fc Hiramatsu| (|2008[); |Matsubara 



(2008) 



(2008) 



Bernardeau et al.| (|2008[ ); jBernardeau fc Valageas 
Anselmi et al.| (|2011[ ) (and more recently Sato & 



Matsubara| ( |2011[ ); |Elia et al.| ( |2011| >; |Wang fc Szalay| ( |2012[ ) 

for the case of biased tracers). Overall the resulting con- 
clusion of these works is that the matter P(k) can be mod- 
eled at the percent level accuracy on weakly nonlinear scales 
(k < 0.2 — 0.4 ft Mpc -1 depending on redshift), improving 
over standard PT. Nonetheless the structure of the solu- 
tions are complex, generally involving a set of couple integro- 
differential equations or multi-dimensional integrals that in 
any case require a time-scale of hours to evaluate. 

In this paper we try to overcome this problem using 
an effective description of the multi-point propagators in- 



troduced in Bernardeau et al. (20081. The multi-point prop 



agators are formally defined as the infinitesimal variation 
of cosmic fields with respect to the initial conditions. For 
Gaussian initial conditions they are equivalent to a measure 
of the cross correlation of final fields with initial configura- 
tions. On large scales, where PT is valid, the propagators 
coincide with the standard kernels in the PT expansion. To- 



wards small scales nonlinear effects drive the propagators 
to zero ( Bernardeau et aL]|2008 l. The full dependence with 
time and scale is then highly non-trivial (besides the fact 
that formally they have a matrix structure). However the 
importance of the MP reside in the fact that they can be 
used as a well behaved expansion basis for equal time cor- 
relators such as the power spectrum, bispectrum, etc. Our 
effective description for the MP accounts only for most grow- 
ing contributions. This, in turn, allows for a rapid evalua- 
tion of the first few terms in the MP expansion of the power 
spectrum. We have thoroughly tested against a large set of 
dedicated N-body simulations both the prescription for the 
different MP (and the matrix structure) as well as the re- 
sulting prediction for the nonlinear power spectrum. For all 
cosmologies investigated we find that our approach is able 
to reproduce N-body measurements at BAO scales at the 
~ 2% level from low to high redshift, with evaluation times 
of about ten seconds. 

This paper is organized as follows. In Sec. [2] we present 
the sets of large N-body simulations of different cosmologi- 
cal models used throughout the paper. In Sec. [3] we briefly 
recall the concept of multi-point propagators and their use 
as expansion basis for the nonlinear matter spectrum. Sec- 
tion [4] goes into more detail with the MP, comparing our 
effective modeling against measurements of two and three 
point propagators in our N-body simulations. In Sec. [5]wc 
use these results to compute the nonlinear P(k) and com- 
pare it to measurements in our fiducial ensemble, with a de- 
tailed discussion of the code performance (evaluation time, 
integration accuracy, etc). Sec. [6] extends this comparison 
to power spectrum measurements at various redshifts in our 
six (6) different cosmological models. Lastly, Sec. [7] contains 
our conclusions. 

We leave for Appendix [X] an important discussion re- 
garding the validity of using PT techniques derived within 
an Einstein - de Sitter cosmology to describe arbitrary 
LCDM models. We carry this out with a novel approach 
using numerical simulations. 



2 N-BODY SIMULATIONS 

We now describe the set of large N-body simulations that 
we developed to test our theoretical predictions. We have 
developed a large ensemble of high statistical significance 
for a fixed cosmological model and we have also carried out 
a set of smaller ensembles for different cosmological models. 
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To this later set we add measurements of P(fc) from some 
of the LasDamas simulations^! 



All simulations used Gadget2 (Springel et al. 20011 



to compute the gravitational evolution, and 2nd order La- 
grangian Perturbation Theory (2LPT) to set up the initial 
conditions at Zi = 49 (Scoccimarro 1998 Crocce et al. 20061. 



2.1 Fiducial Cosmology 

The core testing will be done against measurements in a 
set of 50 N-body simulations, each of comoving box-size 
Lbox = 1280 ft -1 Mpc and 640 3 particles. This set then con- 
stitutes more than 100 h~ 3 Gpc 3 of simulated volume and 
will be particularly important to test our model assumptions 
beyond two-point statistics, i.e. the three-point propagator 
that we discuss below. 

We will refer to this set as the fiducial cosmology (FID). 
The cosmological parameters and relevant information are 
given in Table [I] (see Crocce & Scoccimarro (2008) for more 
details on the simulations). The corresponding particle mass 
was m p = 6x10 h~ Mq . We note that both the mass res- 
olution and the settings employed to run Gadget2 and 2LPT 
ensure that we achieve unbiased measurements of the power 
spectrum at the scales of the baryon acoustic oscillations 



( Heitmann et al.| ( |2010| |) 



2.2 Cosmological Suite 

The main goal of this paper is to provide an efficient pre- 
diction for P(fc) that serves across cosmological parameter 
space. In order to explore this we implemented numerical 
simulations of 5 different cosmological models in addition 
to the FID case, changing the parameters that are of most 
importance to large scale clustering such as matter density 
fl m , spectral tilt n s , linear amplitude of fluctuations as and 
more. Full details are listed in Table [l] 

Each cosmological model was simulated with four (4) 
runs of comoving box-size Lbox = 1250 /i -1 Mpc and 640 3 
particles (Gadget2 and 2LPT settings as in the FID case 
above). 

In addition we use some of the LasDamas simulations 
(McBride et al. 2012, in preparation). LasDamas is a col- 
laborative effort that run 50 boxes of 4 different resolu- 
tions each, tailored to describing the clustering of differ- 
ent galaxy samples in the SDSS-II survey. Here since we 
are interested in large-scale clustering, we only use 4 of 
their largest boxes, named Oriana, each with 1280 3 par- 
ticles within Lbox = 2400 ft -1 Mpc (particle mass m p — 
4.57 x 10 11 ). The outputs for Oriana sample a wide red- 
shift range (z = 0,0.34,0.52,0.97,1.50 and 2.52), allowing 
us to test how the performance of our power spectrum pre- 
scription evolves with redshift. 



2.3 Simulations with Independent Initial 
Positions and Velocities 

Initial conditions in an N-body simulation of a given cosmo- 
logical model are fully specified by initial (random) values 



of particle positions and velocities, i.e. density 5 and veloc- 
ity v perturbation^] Since in linear evolution the vorticity 
component of v decays in time, the only relevant component 
of the velocity field is its divergence 9 = —V ■ v/Hf. More- 
over the linear evolution of 8 and 6 admits two solutions: 
one that grows in time (oc to the linear growth factor) an- 
other that decays away as H(t) in ACDM cosmology. Only 
the first solution survives in the long-time limit. Hence cos- 
mological N-body simulations are always initialized already 
in the "growing mode", that is, setting to zero the decay- 
ing solution by construction. This is achieved in practice by 
requiring that initially (minus) the divergence of the pecu- 
liar velocity field is in phase with density perturbations, i.e. 
<5init ~ #init = <j>- Hence only one scalar field <f> (randomly 
sampled from the linear power spectrum) determines the 
full realization of initial perturbations. 

Here we depart from this standard practice, for the fol- 
lowing reason. Our theoretical framework is built upon the 
multi-point propagators, which can be measured in simu- 
lations as the cross correlations of final and initial fields 
([Crocce fc Scoccimarro l | 2006a I, Bernardeau et al. (20081, 
Bernardeau et al. ( 2012 However if initial fields are in 
the growing mode only certain combinations in the cross 
correlations can be tested. For instance, the 2-point cross- 
correlation r a b = (ab) (related to the 2-point or nonlinear 
propagator G a b) should a priori have 4 independent combi- 
nations out of a — (Sf,0f) and b — (di,0i); but in practice 
only the projection along Si = 8i is measurable if the ini- 
tial conditions are in the growing mode: G a = G a bUb = 
Gaii +G a 8i, where Ub = (1, 1) indicates that the initial con- 
ditions are given by (5i,0i) = (1,1)0. Thus, to date only 
the "density" Gs or "velocity" Go two-point propagators 



have been tested against simulations ( Crocce & Scoccimarro 
2006a| [Bernardeau et al.|2008| [2012| 



In order to measure all the components of the propaga- 
tor, and hence test our analytic predictions in much more 
detail, we have performed for the first time simulations with 
independent particle positions and velocities, leading to in- 
dependent density and velocity perturbations. This was done 
as follows. We first generate two Gaussian random fields, 
<f>\ (k) and <j>2 (k), out of the same linear power spectrum Po. 
We then run two simulations mixing the initial conditions 
in densities and velocities: 



run 1 



&nit(k) = . 

ftnit(k) = 



'i(k) 
> 2 (k) 



(1) 



http:/ /lss. phy.vanderbilt.edu/lasdamas 



and the opposite with run 2. In this way each run has initial 
conditions which are general in their initial values of density 
and velocity perturbations (hence a linear combination of 
growing and decaying mode solutions), such that Ps^ (fc) = 
P Mi (fc) = Po(fc) while P M ,(fc) = 0. 

At the practical level the procedure described above is 
equivalent to use any standard initial condition generator to 
produce two random sets of initial positions and velocities 
(out of the same linear Po), and then combine the initial 



7 By velocity perturbations we mean the velocity field after the 
Hubble flow at the given position has been subtracted off (the 
so-called peculiar velocities). 

8 We are assuming Gaussian initial conditions in this pa- 
per, otherwise the relationship between propagators and cross- 



correlations is more complicated, see Bernardeau et al. (20101 
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Figure 1. Two-point (nonlinear) propagator for the density field: model vs. measurements in N-body simulations at z = 0, 0.5, 1. The 
model in Eq. | |16| l performs remarkably well at all redshifts shown. The dashed line shows the corresponding high-fc limit (which is only 
reached at very high-fc, not shown here). Lower panels show the ratio of the measurements to the two different analytic descriptions and 
stress that the accuracy of Eq. JTm is at the percent level. 



positions of the first set with the initial velocities of the 
second to perform the first simulation run (and viceversa 
for the second run). See for instance Scoccimarro (19981 for 



a more detailed explanation on how initial positions and 
velocities are set in simulations. 

To some extent these two runs are not fully independent 
from each other, since initial positions and velocities are 
interchanged. We repeat the process twice so we end up 
with four runs of simulations with such initial conditions, 
our measurements below average over such four realizations. 



3 MULTI-POINT PROPAGATOR EXPANSION 

In this paper we will work with the so called multi-point 
propagator expansion of (equal-time) correlators of cosmic 
fields (Bernardeau et al. 20081. In particular we are inter- 



ested in the nonlinear density power spectrum, which in this 
framework is given by 



P SS (k,z) = J2 rl / M k -*h...r) [r< r) (q 1) ...,q r ;z; 



X P { qi ) . . . P (q r ) d 3 qi ...d 3 q r 



(2) 



Here Po denotes the spectrum of perturbations at some 
initial time and are the multi-point propagators, de- 
fined as the (ensemble averaged) variation of late time cos- 
mic fields ^ a = (5, 6) with respect to the initial conditions 

(f> a = ^a(zi), 



1 , (Tfrafos) 
7-1 \ 



i>i( k i 



(k, 



foCk-kr. 



(ki 



, k r , z) 



(3) 



where ki... r = ki + . . . + k r . In the most general scenario 
Eq. Q should allow for arbitrary initial spe ctra of density 



and velocity fields (|Bernardeau et al. 2008 J. But for 



simplicity we have assumed that initial fields are adiabatic 
and in the growing mode, thus P„& = u a utPo(k) with u a = 
(1,1). This translates into the more compact expression for 
the propagators 



= r 



(r) 



«6i 



(4) 



used in Eq. |2| evaluated for the density fluctuations (a = 1, 
or a — 5). Equivalent expressions to Eq. (J2J hold for Pee 
(replacing Ts by Te) and Pse (using a cross-term r^re). 

Equation |2| results from the resummation of a whole 
set of (infinite) terms in the standard PT expansion of P(k). 
Unlike the standard approach, it is a sum of positive terms 
each of which dominates only in a narrow range of scales. 
Notice that now each multi-point propagator has contribu- 
tions to all orders in PT, and depend also on Po- At low k 
they can be described perturbatively while their asymptotic 
properties at large k can be computed beyond perturbative 
expansions. However, to implement in practice Eq. |2| one 
must have a description of multi-point propagators at all 
scales (and times), matching the perturbative calculations 
at low k to the resummed asymptotic behavior at high k. 
Achieving this in a way that is fast and accurate enough for 
the needs of cosmological surveys is the goal of our work. 

The numerical evaluation of increasing orders in Eq. |2| 
becomes very demanding rather quickly. In order to main- 
tain a fast evaluation time we will concentrate on quasilin- 
ear scales and implement Eq. |2| up to r = 3 for which we 
need a description of the two, three and four-point propaga- 
tors. In what follows we focus in discussing the multi-point 
propagators in more detail, with emphasis on our particular 
description of r (1) ,P (2) and r (3) . We then discuss in Sec. ^ 
the computation of the power spectrum and comparison of 
our predictions against measurements in simulations. 
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4 THE MULTI-POINT PROPAGATORS 

The two-point or nonlinear propagator, first introduced by 



Crocce & Scoccimarro (2006b) in the context of Renormal 



ized Perturbation Theory, is the (ensemble averaged) vari- 
ation of late time cosmic fields ^ a = (<$, 9) with respect to 
the initial conditions <f> a = ^ a (zi), 



<5 D (k - q) Gab [k,a) = ( . . ), 



b(q) 



where we adopted the notation G a b 



and used as time 



variable the growth factor 

gQ This object emerged from the 
resummation of an infinite subset of contributions to the per- 
turbative expansion of the power spectrum. This effectively 
"renormalized" the linear growth factor into a fully nonlin- 
ear and scale-dependent function: the nonlinear propagator 
Gat- The precise way in which it is a renormalized version of 
the growth factor can be seen for Gaussian initial conditions, 
in which case G a b fully describes the cross-correlation be- 
tween initial and final conditions, i.e. (^a^b) = G ac (<^ c <^t,}. 
Also, one can easily show from Eq. Q that an expansion 
of cosmic fields in terms of their initial values (as done in 
standard PT) leads to 



G a b = g a b + "nonlinear (loop) corrections" 



where g a b is the standard linear propagator, 



-3/2 



(6) 



(7) 



Hence, on large scales were linear PT becomes a good ap- 
proximation we recover 



GabUb — > a, as fc 







(8) 



while on small scales nonlinear effects become dominant 
driving G to zero, as initial (linear) and final fields become 
decorrelated. In Crocce & Scoccimarro (2006a I it was shown 



that in this limit it is possible to resum all the dominant per- 
turbative orders exactly, leading to 

Gab(k, a) « g a b(a) exp[-^k 2 al] as fcer d > 1, (9) 

where the characteristic scale of decay is given by the r.m.s. 
one-point displacement field that to most growing order co- 
incides with the amplitude of large-scale velocity flows, 



2 



{a -If fP . 



-d q. 



(10) 



In Bernardeau et al. (2012 1 it has been explicitly shown 



that this result can be derived in a very general framework, 
the eikonal approximation, irrespectively of the time depen- 
dence of the large-scale flows. This partially extends the 



result of Anselmi et al. (20111 who were able to resum a 



9 In what follows we will assume the structure of the theory is 
that of an Einstein de Sitter universe (Q m = 1) for which the 
growth and scale factor coincide. The validity of the calculations 
is not significantly affected on large scales by this assumption if 
we replace a by the appropriate growth factor D+ for the cor- 
responding ACDM model. We test and discuss this in detail in 
Appendix [A] This approach is equivalent to an approximation 



about the linearly decaying modes, i.e. D- 



D 



-3/2 



sub-leading set of perturbative contributions, what led to a 
slight modification of the damping scale Od . 

The importance of the high-fc limit asymptotics in 
Eq. |9| stands from the fact that in RPT all diagrams 
for correlators now have the linear propagator replaced 
by the renormalized propagator G a b inside all loops. As a 
result, after this resummation, large scale predictions are 
least sensitive to what is going on at small highly-nonlinear 
scales, where even the pressure-less perfect fluid approxima- 
tion (used to derive these results) breaks down. This is in 
contrast with standard PT where, depending on the num- 
ber of loops, the propagator has the wrong asymptotics, 
G a b(k — > oo) — > ±oo. 

The two-point propagator (sometimes also called the re- 
sponse function) turned out important not only to RPT but 
also in other resummation schemes such as the Path Integral 



approach of Matarrese & Pietroni (20071, Closure Theory 
( Taruya fc Hiramatsu | |2008| ) , large-N expansion (|Valageas 



cently the extended TimeRG of Anselmi & Pietroni (20121 



2007a I, Lagrangian schemes ( jMatsubar a 2008) and more re 



Note, however, that in some of these cases the propagator 
behavior at high-fc can be different than in RPT and what 
we present here (which are in good agreement with simu- 
lations). In the Closure and Large-N cases, there are un- 
physical oscillations superposed with decay, whereas in the 
Lagrangian approach of Matsubara (2008 1 part of the prop- 
agator remains perturbative, and thus the same issue about 
violation of high-fc asymptotic as in standard PT occurs. 
Another key step forward was achieved in |Bernardeau| 



et al. ( 2008 \ were it was shown that the concept and re- 



sults of the two-point propagator could be extended to an 
arbitrary number of points, 



(5 p *„(k,a) 



r> = 



p\ x <5-/> bl (ki)...(50 i , p (kp)' 

fo(k-k x ..,) r&> ..^(kr,...,^), (11) 

where ki... p = ki + . . . + kj^j In this case, a relation anal- 
ogous to Eq. (|6| is obtained 



ri™' = J- a n> + loop corrections 



(12) 



-O) 



where on large scales one recovers the well known 
a n (F n ,G„) kernels in PT (assuming growing mode initial 
conditions and keeping only the fastest growing contribu- 
tion), 

Pi n) ~a n {F„(k 1 ,...,k n ),G„(k 1 ,...,k„)} as fc^O 

for a = 1, 2 (density or velocity divergence fields respec- 
tively). Again, on smaller scales ri™' is expected to be 
driven to zero, as the p-point propagator can be shown to 
be proportional for Gaussian initial conditions to the cross- 
correlation (\£ r <^>. . . 4>)/Pq , a generalization of the two-point 
result. Remarkably, the multi-point propagators also admit 
the resummation of the dominant behavior in the high-fc 
limit, yielding, 



T ( a n) exp[- 



1,2 2i 
2 fe a d\ 



(13) 



10 Note that r( p ' denotes the (p + l)-point propagator, by trans- 
lation invariance it depends only on p wavenumbers in Fourier 
space. 
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in strict analogy to Eq. jjjl. Because of this, Eq. |2| for 
the power spectrum can be thought of as an expansion in 
terms of cross-correlations that are i) always positive and 
each term dominates in a narrow range of scales, ii) are well- 
behaved even in the nonlinear regime due to Eq. (13 1. The 



convergence of such an expansion in the nonlinear regime has 
been verified for the case of the Zel'dovich approximation, 
where the exact result is known and the expansion can be 
carried out to a large number of loops (JValagcas 2007b| ). 

So far we have discussed limiting expressions of the 
multi-point propagators in the low and high-fc regimes. To 
continue we need prescriptions valid at all scales that can 
be integrated over as in Eq. This will be the subject of 
the following sections. 



4.1 Two-point (nonlinear) Propagator 

We now discuss our prescription for matching the small and 
large fc regimes of G a b and thus reconstruct the full two- 
point propagator. The procedure is similar to that in RPT 
(Crocce & Scoccimarro 2006a) but is simplified because we 
are interested in describing the propagator at late times (un- 
like RPT, we will not be doing time integrations over the 
propagator here). 

The next-to-leading order corrections to the linear prop- 
agator are given by a one-loop computation. If we neglect all 



sub-leading time dependencies the result reads (see Crocce 



& Scoccimarro (2006a I for the full expression otherwise), 



5G lloop = 


^a 3 f(k) + 0(a 2 ) 




8G{ 2 oop = 


2 -a A f{k) + 0{a 2 ) 




<5G*2i oop = 


^a 3 g(k) + 0(a 2 ) 




<5G22°° P = 


\a A g{k) + 0{a 2 ) 




(14) 



where 
f{k) 



1 



504fc 3 g 5 



\6k 7 q - 79k 5 q 3 + 50<? 5 fc 3 - 21fcg 7 



+ 



-(fc 2 -g 2 ) 3 (2fc 2 + 7g 2 )ln 



k + q\ 2 



+ 



^ [6*^-41^ + 2* 

3(fc 2 - g 2 f(2fc 2 + g 2 )ln^ 



Po(q) d 3 q, 
f - 3kq 7 
l]Po(q) d 3 q. 



Notice that both a 2 f(k),a 2 g(k) -> -k 2 aj/2 for ka d > 1. 
Therefore we can put together Eqs. ( 14 1 and the most- 
growing term in Eq. Q and "exponentiate" the propagator 

as, 



Sll + 5Gii oop - 


+ Gn = 


3 2 

-oexp[o f(k)} 




gi2 + 8G^ op - 


+ Gl2 = 


jjaexp[a 2 /(fc)] 




921 + 8G 21 p - 


+ G 2 i = 


|aexp[a 2 g(fe)] 




322 + SG 22 p - 


-> G22 = 


jjaexp[a 2 #(fc)] 


(15) 



quasi-linear scales as well as the dominant large-fc asymp- 
totic (at late times). They do not exactly match the ones 
given by Crocce & Scoccimarro ( 2006a[ ) for RPT. Those ex- 
plicitely preserved sub-leading time dependencies needed to 
correctly integrate the evolution from initial conditions. 

For initial conditions in the growing mode, cj> a oc u a — 
(1, 1), the relevant quantity is the two-point "density" prop- 
agator Gs ~ GiaUa = Gn + G12 that simply reads, 



G s {k,z) = D+(z)e X p [f(k)D 2 + (z)] 



(16) 



In Fig. [T] we show how this prescription performs against 
measurements of the propagator in our FID ensemble of N- 
body simulations (top entry in Table if"" Left, middle and 
right panels correspond to redshifts z = 0,0.5 and 1 respec- 



tively, the model from Eq. ( 16 1 is shown in solid blue while 
the high-fc limit expression in Eq. |9| in dashed line. De- 
picted error bars correspond to the variance over the en- 
semble. This prescription, although simple, describes the 
measurement at the sub-percent level at all redshifts and 
for all scales of interest (bottom pa nels in Fig. [T] show the 
corresponding fractional errors). In Taruya et al. (20121 a 



thorough analysis of the performance PT predictions at one 
and two-loop order is performed. It is found that for z> 1 
two-loop results can improve upon one-loop predictions. We 



stick however to one-loop prediction as given in Eq. ( 16 \ in 
this work as it gives predictions of ample precision in all 
regimes of interest. As we will see in Sec. [6] this holds true 
not only for our fiducial (FID) cosmology but for all the 
cosmological models listed in Table [l] 



4.2 The Full Matrix Structure of the Propagator 



As discussed in Sec. |2.3| in standard cosmological simula- 
tions initialized in the growing mode it is only possible to 
measure two linear combinations of the two-point propaga- 
tor. In order to probe the full matrix structure (the four 
elements of G a b) we must initialize simulations with two in- 
dependent random Gaussian fields for particle positions and 
velocities, respectively. Let us now discuss how the resum- 
mation in the high-fc limit is changed when we have inde- 
pendent initial density and velocity fluctuations. 

A good starting point is the analysis of the structure of 
the one-loop contribution to G a b calculated from, 
rv rsi 

SG ab° V ~ 4 / dsi / d,S2g ac (r)- si)7 de(k,q, k- q) 



jo jo 

x gdf(si)g eg (si -s 2 ) x 7 sM (k + q',q',k) 

x 9hj(s2) <0/(q)</>j(q')) 9ib(s 2 ) (17) 

where 7 a (, c is the vertex function in standard PT (see for 
instance Eq. (22) in Crocce & Scoccimarro (2006a I for a 
general derivation of this expression), since this is what then 
gets exponentiated by the resummation procedure. 

In the standard growing mode case the correlator of 
initial conditions is, 



>/(q)0i(q')) = w/WjPo(q)&(q+q') 



(18) 



11 The measurements have been performed by cross-correlating 
the density field at the desired redshift with itself at the initial 



These expressions recover the correct one-loop result at 



conditions (set in the growing mode), see |Crocce &: Scoccimarro| 
1 2006a 1 for a detailed discussion about the estimator. 
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Figure 2 . Components of the Nonlinear propagator. We show for the first time the four individual components of the nonlinear propagator 
(normalized to the linear growth factor) measured in dedicated simulations with independent <5 and 9 initial conditions. In our model the 
decay of the density propagators G gg and Ggg is given by exp[(13/25)/(fc)-D 2 (,z)] while the velocity components Ggs and Ggg are given 
by exp[(13/25)g(fc)-D 2 (,z)], see Eqs. 1 15|22 I. Dashed lines show for reference the decay obtained in the high-fc limit exp(-fc 2 (13/25)o"^/2) 
(same for all). 



with u = (1,1). For our mixed mode initial conditions we 
have instead (see Sec. 2.3 1, 

( ^(q)^(q') ) = <5£P (q)fo(q + q'), (19) 

where the Kronecker symbol 5fj = 1 if / = j and other- 
wise. This means that instead of evaluating 



94f(si)g h j(s 2 )ufUj = e° 1+ " 2 u d u h 



(20) 



in the standard case (with s — In a), we have to compute 
13 



L + 32 



25 



UdUh + decaying mode, (21) 



-3(si+s 2 )/2 



where the decaying mode piece evolves as e 
Therefore, neglecting the decaying mode contribution we see 
that the overall effect of using independent random field ini- 
tial conditions is to renormalize 

13 



{/(*).»(*)} 



25 



{/(*),«/(*)} 



(22) 



in Eqs. (14 1. This can be carried out to all orders leading to 



the same high-fc limit resummation as Eq. ( 16 \ except that 

(23) 



2 
0~d 



13 2 
25^ 



Thus our model for the full propagator for independent S 
and 9 initial conditions is simply the one in Eqs. ( |15| ) with 
the replacement given in Eq. ( 22 l . 



In Fig. [2] we show the four different components of the 
propagator measured in the simulations with independent 
S and 9 initial conditions against predictions by the above 
model (in solid blue) and the high-fc asymptotic (in dashed 
black). The propagator was measured following, 

O Mk)fo(-k)) 

Po(fc) 

where now the initial conditions <f>b are different for density 
(b — 1) or velocities (b = 2) according to Eq. |TJ. The veloc- 



G ab (k) = 



(24) 



ity divergence fields in Eq. ( 24 1 were estimated following the 
procedure describe in Scoccimarro (20041. Reassuringly all 



the four components follow the expected theoretical decay 
towards small scales. This is an interesting and important 
cross-check, in particular for resummation schemes such as 
RPT or Closure Theory, that integrate the individual com- 
ponents separately rather than the "density" and "velocity" 



propagators that until this paper were the only combinations 
tested against simulations. 



4.3 Three and Four-Point Propagators 

The three-point propagator was introduced by Bernardcau 



et al. 



et al. 



( 2008 1 and recently studied in detail in Bernardeau 



(20121. In the later work a general scheme (called 



RegPT) to interpolate between small and large scale infor- 
mation for any propagator was proposed. In particular, the 
RegPT prescription for the three-point propagator is given 

by, 



n(2)RcgPT 



(kl,k 2 )= [i 



(2)Trcc 



(ki,k 2 ) + $r 



(2)onc — loop 



(ki,k 2 ) 



1 1 , 2 2 r ,(2)Trcc/ 1 1 \ 

+ 2 fc adT ( k i» k 2j 



cxp(-fc 2 <Td/2) 

(25) 



where k = ki + k 2 . Here r (2)Trcc , 5r (2)onc - loop and a d de- 
pend on time. The one-loop term in this expression is de- 



scribed in detail in Bernardeau et al. ( 2012 1 and involves one 
integral over Po(q) for each triangle configuration (ki,ka). 
Although Eq. |25| gives a very good agreement with mea- 



surements in simulations we have found that its usefulness to 
compute the one-loop P(fc) is limited because it takes a long 
time to evaluate. We hence seek an alternative prescription. 



We have found that, in analogy to Eq. ( 16 1, the follow- 
ing expression (k = ki+k2) 



rf>(ki,k 2 ;z) = D^(^)P 2 (k 1 ,k 2 ) eMf(k)D'+(z)} (26) 



yields very similar results to that in Eq. ( 25 \ and virtually 



the same one-loop power spectrum after the correspond- 
ing momentum integration. Figure [3] shows measurements 
of r' 2 ^ for different triangle configurations together with 
the prediction by the model in Eq.(|26[) in solid black (used 
throughout this paper) and RegPT from Eq. ( 25 I in dashed 
black. Left and Middle panels correspond to equilateral with 
fci = fc2 = fc3 = fc (top left); colinearwiih k\ = fc2 = fc/2 and 
fc3 = fc (top center); elongated (bottom left) with fci = fc/2 
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Figure 3. Nonlinear three-point propagator: analytic predictions vs. measurements. Different panels show measurements of T( 2 ' (qi, <?2j 93) 
in our fiducial ensemble of simulations for different triangular configurations, as indicated in the y-axis label (see text for details). Solid 
line corresponds to the interpolation scheme proposed in this paper, see Eq. J26J. Dashed line to the one introduced in Bernardeau 
et al. (2011) (RegPT), see Eq. l |25[ l. They are mostly indistinguishable for almost all configurations and agree with the measurements 
remarkably well. In particular for right panels where we show the configuration that contribute the most to the one-loop power spectrum 
at (fixed) wavenumber k = 0.06 ft -1 Mpc (0.1 /i -1 Mpc). Error bars correspond to the variance over the ensemble and results (in left and 
middle panels) are plotted against k% = k. 



and &2 = &3 = k; and squezeed (bottom center) fci = fe/4 
= kz — k configuration j^j 
Notice that the theory is binned in the same way as 
the data, this is essential to recover the correct asymptotic 



behavior at low-fc (see Bernardeau et al. (20121 for details 
on the r' 2 ' estimator and the binning correction). 

From these panels it is clear that both models perform 
very well for all these configurations with a slight over- 



prediction by Eq. ( 26 1 for squeezed configurations. 



In addition the right panels of Fig. Q show the same 
comparison for the configuration that would yield the dom- 
inant contribution to the one-loop computation of P(k). 
From Eq. (2| we see that the one-loop power spectrum is 



12 We assume the final (nonlinear) density field has wave-vector 
&3 (last argument). Hence the three-point propagator is only sym- 
metric with respect to the 1st and 2nd indices that corresponds 
to the initial (linear) fields. 



of the form 



5 Hoop 



(fc) 



4tt 



Po(qi)qidqi / ^0(92)^2^52 



x [rf (91,92,*!)]' 



(27) 



Hence by symmetry reasons the most relevant configura- 
tion for a given k is roughly T^(q, q, k). Figure ([3J shows 
this configuration for k — 0.06 h Mpc -1 (top right panel) 
and k — 0.1 h Mpc -1 (bottom right panel). Here again the 



model in Eq.(26l describes the N-body results remarkably 



well yielding the same answer as RegPT (notice the dashed 
and solid line are on top of each other). 

The four-point propagator is basically a measure of the 
trispectrum between final and initial density fields. Thus it is 
difficult to measure from the N-body simulations. Nonethe- 
less from theoretical grounds we do know the behavior at low 
and high-fc, and have no reason to expect a different behav- 
ior at intermediate scales from the one already probed for 
the two and three-point propagators. Hence we will adopt 
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the following prescription (k = ki.23), 

If>(ki,k3,k8;*) = D^zJFstki.ka.ks) exp[f(k)D 2 + (z)], 

(28) 



for the four-point propagator, in full analogy to Eqs. ( 16 \ 
and ( |26[ ). This prescription will then satisfy the low-fc and 
high-fc asymptotics. 

Provided with prescriptions for the propagators up to 
four points, we are now ready to discuss the multi-point 
expansion for the power spectrum. 



5 POWER SPECTRUM 

In this section we present the prescriptions we adopt to 
do the actual computations of the power spectra. We call 
MPTbreeze this implementation and comment at the end of 
the section on possible alternative approaches. 

To describe the power spectrum at mildly nonlinear 
scales we implemented the first three terms in the expan- 
sion in Eq. Q. In diagrammatic language they correspond 
to renormalized versions of tree level, one and two loops 
respectively (by renormalized we mean that the renormal- 
ized propagators include themselves loops to all orders). The 
tree-level term is simply given by 



P trcc = [r (1) (M)]"P 



(29) 



and c oincides with the propagator r enormalization term of 



RPT (Crocce & Scoccimarro 



2006bh. Here T^ 1 ' is given by 



Eq. qisp and Po is the linear, post-recombination, spectrum 
of fluctuations. This term is the one that contain most in- 
formation on narrow band features of the primordial per- 
turbations such as the Baryon Acoustic Oscillations. The 
"mode-coupling" contributions start with the next term in 
Eq. |2} that reads, 

J\-loop(*, z) = 2j d 3 q[rf } (k - q, q; z)] 2 P (|k - q|)P ( g ) 

(30) 

Assuming k along the z-axis it can be easily turn into 
Pr-joopCfc, z) = 4arf dx J dq[Vf (p, q, y; z)] 2 P (p)P (q) 

where x = k • q/(fc q) and we have introduced p = k — q and 
y = p • q/(pq) for clearness. Here is given by Eq. (261 

with F 2 (p, q, y) = 1 + 1(^ + 1,) + ^ the standard 2nd-order 
PT kernel (e.g. |Bernardeau et al.| (|2002[)). We are then left 



with an integration in 2-dimensions which can be rapidly 
evaluated with a standard Gaussian quadratures routine. 

The third term in Eq. |2| is slightly more involved but 
can be treated similarly, we now have 

P2-ioo P (M) = d 3 ^ J d 3 q. 2 [rf ) (\t-q 12 ,q 1 ,q 3 -,z)] 2 

Po(|k-q 12 |)P ( 9l )Po(g2). (31) 

Here we can take k along the z-axis and ^ in the x — z, 

k = k (0,0,1) 

q± = gi(sin#, 0, cos 9) 

•32 ~ 92(sin <£sin a, cos <f>sina, cos a), 




0.05 0.10 0.15 0.20 0.25 
k [h MpcT 1 ] 

Figure 4. The multi-point propagator expansion presented in 
this paper (solid blue line) against measurements of P(k) in our 
FID ensemble of N-body simulations (top entry in TableQ at z = 
0, 0.5 and 1. The dotted red line is linear theory and solid black is 
halof it. The evaluation time of the multi-point expansion shown 
in each panel is at most five seconds. 



so we are then left with the following 5-dimensional integral 
P2-ioo P (fc, z) = ^'"■j ^l 1 J d 1 2 j dx j dy j d(j> 

[I? ) (q 3 , Qi , q 2 ; *)] 2 Po (9s)Po (gi)Po (32) 

were we introduced q 3 = k — q 12 ; x = cos 6 = k- <li/(kqi) 
and y = cos a = k-q 2 /(fcg2) are integrated in [—1, 1] and <j>, 
the azimuthal orientation of q 2 , in [0, 2ir]. The propagator 
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rf> is given by Eq. |28| in terms of f(k), Eq. |15j), and F s , 
which in turn is solved iteratively in terms of (F2, G2) (see 
" 20021"). 



Bernardeau et al. 



e.g. 

It should be noted that all the required integrals safely 
converge for any realistic shape of the matter power spec- 
trum. The regularization scale <Td is finite as soon as the 
power spectrum is steeper than fc _1 in the low wave-mode 
limit (n > —1). All terms involved in the computation of 
the power spectrum actually require the same conditions 
in the IR domain. The convergence in the UV domain is 
more diagram dependent. The condition for the existence 
of <7d demands that n < — 1 if n is the power spectrum 
index at large wave-modes. The convergence properties of 
fi-ioo P (fc, z) and P2-ioo P (fc, z) are all determined by the be- 
havior of the 7a6c(k, q, k — q) (symmetrized) vertex func- 



tions, introduced in Eq. (17 1. They obey the scaling relation, 



7 a6c (k,q,k-q) 



(33) 



when 5 > fc. As in the expressions of Pi_i oop (/c, z) and 
P2-ioop(fe, z), two such factors are introduced in the high 
q limit, the convergence of the former is secured as soon as 
n < 1/2 and that of the latter when n < —2/3. It therefore 
does not lead to new constraints besides the existence of aa- 
This is at variance with the convergence properties of the 
multi-loop corrections to the two-point propagators. They 
indeed lead to much stringent constraints in the UV domain 
as discussed in |Taruya et al.| ( |2012[ ). We will comment fur- 
ther on the consequences of those properties in Sec. |5.2| 
To perform the integration of Eq. ( 32 I we employ a 



MonteCarlo algorithm routine called Vegas (Lepage 1978 
|1980[ ), included within the vl.4 CUBA library of multidi- 
mensional numerical integration routines describ ed in|Hahn| 



( 2005[ f 



5.1 



This is discussed in more detail in Sec 
We are now ready to compute the model prediction 
and compare it with measurements in N-body simulations. 
Figure [4] shows the two-loops multipoint expansion in solid 
blue line, linear theory in dotted red and halofit (Smith 
et al.|[2003[ ) in solid black, at z = 0,0.5 and 1 (top to 
bottom panels). Symbols with error bars are the corre- 
sponding measurements of P(k) in our fiducial cosmologi- 
cal model (error-bars correspond to the variance over the 
ensemble of 50 simulations, see Table [y. All lines are di- 
vided by a smooth broad-band linear spectrum to reduce 
the y-axis dynamic range. Overall, the multi-point expansion 
match the measurements at the 2% — 3% level up to scales 
k ~ 0.16/iMpc- 1 ,0.19/iMpc- 1 ,0.23/iMpc" 1 at z = 0,0.5 
and 1 respectively. These values coincide roughly with 1 
at the given redshift, with ad given by Eq. (10 1. Notice also 



that these scales are slightly beyond the Baryon Acoustic 
Oscillations region. Improving upon this ~ 2% accuracy re- 
quires a more precise ansatz for the multi-point propagators, 
in particular for F*- 3 ' because it dominates within this region 
(k > 0.1/iMpc" 1 ). 

We have also investigated whether the addition of the 
next term in the expansion of Eq. |2| , that corresponds to a 
3-loop computation, extends the agreement to higher wave- 
modes. The result is that it does but only slightly because 



of the already strong exponential suppression of the corre- 
sponding five-point propagator. In turn, the numerical eval- 
uation of the 3-loop integral becomes very lengthy. An alter- 
native path is to combine our PT approach with halo model 



prescriptions, in the spirit of Valageas & Nishimichi (20111 



This will be the subject of future work. 

For our FID cosmology we have checked that the recent 



ansatz by Tassev & Zaldarriaga (2012) under-estimates the 
N-body measurements by 4% at z = 0, gives a ~ 1% match 
at z = 0.5 and over-estimates P(k) by ~ 2% at z = 1 
(at BAO scales). This ansatz is based upon a split between 
long and short wave- modes, which was somewhat arbitrarily 
choosen to be a sharp fc-cutoff at A = k/2. The strength 
of nonlinear corrections in the model, hence also the above 
residuals, depends (systematically with redshift) on A. 



5.1 MPTbreeze: Code performance and convergence 

The major advantage of the method presented in the pre- 
vious section is that the evaluation time is of the order of 
a few seconds, comparable to that of linear Boltzman codes 
used to compute the transfer function of different species. 
This is opposed to resummation techniques such as RPT 
and Closure Theory that take significantly longer times to 
compute. Hence our approach it is very well suited for sam- 
pling the large-scale structure likelihood of present and fu- 
ture datasets at BAO scales. 

In addition note that due to the structure of the expan- 
sion the same evaluation done for a given redshift can be 
properly re-scaled to another one at almost no extra cost, 
just recomputing G and the corresponding growth factor. 

It seems then appropriate to discuss the performance 
of the code in more detail and the numerical convergence. 
Clearly, the component that is numerically intensive to eval- 
uate is the 5 dimensional "2-loop" integration given in 



Eq. (32 1. As mentioned before we have perform it with the 



Monte Carlo algorithm Vegas that uses importance sampling 
as a variance-reduction technique. 

The accuracy of the integration, and hence the evalua- 
tion time, is controlle d by setting the required absolute e a b s 
and relative e re i error if^( |Hahn||2005| ). 

For s re i — 1% (with e a bs = 0) the full computation 
of the z — 1 power spectrum including tree, one and two 
loops in steps of 5k — 0.005 h Mpc" 1 (roughly the fun- 
damental mode in our simulations) up to k ma x = o"^ 1 = 
0.244 h Mpc -1 (49 bins) evaluates in only 5 sees. The same 
integration with e re ; = 0.5% requires 10 sees and roughly 3 
minutes for e re i = 0.1%. Moreover, we find that generally 
the three accuracies yield the same power spectrum to a sub- 
percent level thus e rc i = 1% (the fastest) seems a reasonable 
choice. 

A similar test at z — 0.5 is even quicker because the 
validity of the expansion is more limited. Now k max — °"d 1 — 
0.195 ft Mpc -1 (39 bins) evaluates in just ~ 3 sees, setting 
£rci = 1% (7 sees, with e re i — 0.5%). Notice that all the 
timings reported are for single CPU (and the code compiled 
with the Intel compiler ifort). 

In all cases the numerical integration in momentum 



Publicly available at http://www.feynarts.de/cuba/ 



14 Convergence is achieved once the estimation I of the integral 
I satisfies \I — I\ ^ max(e a i, s , e IC \I) 
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Figure 5. Nonlinear two-point propagator at z = for different cosmological models. In each panel symbols with error bars correspond 
to the measurements of the propagator over four simulations of the given cosmology. In solid red lines we show the prescription used 
throughout this paper corresponding to the exponentiation of the most-growing one-loop contribution: exp [D-^(z) 2 f(k)~\. It agrees 
with the measurements at the sub-percent level for all the scales of interest in all cases studied. 



space is done from a fixed q m m ~ 10 -4 up to some cut-off 
scale q c . A priori the box-sizes of our simulations are large 
enough that wave-modes longer than L^ ox have negligible 
amplitude to alter the measured power spectrum in the sim- 
ulations. Hence we assume no finite box-size effect and take 
q-min arbitrarily small. On the other hand we find that we 
need q c ~ 1/iMpc -1 in order for the integrals to converge 
within 0.5%. This mild sensitivity to the ultraviolet (UV) 
cut-off is due to the particular prescription for MP adopted 
here, with the decay standing as an overall multiplicative 
function (T^ n > = G(k) x F n ) that factors out of loop inte- 
grals such as Eqs. ( 30|31 I. This is unlike approaches such 
as RPT where the full nonlinear propagator is integrated 
inside all loops, what screens much strongly the UV regime 
where e.g. shell-crossing enters (and convergence is achieved 
already for q c ~ 0.5 h Mpc -1 ). 

Lastly we have also performed our integrals with other 
Monte Carlo algorithms implemented within the CUBA li- 
brary, reaching always the same answer obtained with Vegas 
but employing more time to converge. For example Suave 
(importance sampling with globally adaptive subdivisions) 
employs 1 minute 25 seconds to evaluate the 2 = 1 spectrum 
for e rc i = 0.005. Divonne (Friedman & Wright 1981a|b l, 
which uses stratified sampling, requires 61 seconds for the 
same task. As described above, Vegas is substantially faster 
requiring only 10 seconds (or less for higher £ rc i). 



5.2 MPTbreeze compared to RegPT 

It is possible to construct alternative implementations of MP 
resummations. RegPT is an alternative proposition that relies 
on slightly different choices for the computation of the multi- 



point functions and on numerical implementations (Taruya 
|et al.|[2012[ ). Let us first stress that to a large extent both 
approaches share the same advantages and disadvantages: 
the power spectrum is constructed out of a sum of positive 
terms and each of these terms could be computed on its own; 
the resulting power spectrum exhibits a large-fc cutoff which 
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Figure 6. Same as Fig.\5\but at z = 1 (and for three cosmological 
models). The simple model for Gs (solid red line) still performs 
at the sub-percent level in all scales and models. 



signals the limit of the validity range of the computation. 
The main difference in the implementations comes from the 
fact that, in the RegPT case, all two loop order terms are 
taken into account following the prescription proposed in 
|Bernardeau et al~| ([2012 ) in the computation of the prop- 
agators. This is not the case in this work where the loop 
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Figure 7. Power spectrum performance for different cosmological models at z=0 . Top left panel shows the linear power spectra in each 
of the simulated cosmologies divided by a smooth broad-band power. Remaining panels show measurement of P(k) over four independent 
simulations for each model. We used the same P smoo th (except for its normalization) for all panels. In these panels dashed red lines 
correspond to linear theory, black solid lines to halof it and solid blue lines to the the multi-point expansion presented in this work (see 
text for more details). The later is found accurate at the ~ 2% level on BAO scales throughout!! all cosmologies investigated. 



corrections to the propagators are computed with a more 
phenomenological approach. 

Because of these differences, the computational difficul- 
ties of th e tw o approaches vary. The computation of the ex- 
pression (301 when the one-loop correction to F' 2 ' is taken 



into account reveals quite costly for a direct computation. 
The MPTbreeze implementation avoids this difficulty. Even 
with direct Monte-Carlo integrations it is then possible to 
obtain the results within a very short timj^] 

The other main difference concerns the converging prop- 
erties of the involved diagrams. As shown in Sect. [5] dia- 
grams involved in the MPTbreeze computations have good 
converging properties. In practice, for wave modes of the or- 
der of 0.1—0.5 ft/Mpc the results depend on wave modes that 
are of comparable scales. This is not necessarily the case for 
the RegPT implementation. As it will be stressed in a forth- 
coming paper , the two-loop corrections to the two-point 
propagator are sensitive in particular to wave-modes well 



15 In |Taruy a et al. ] l |2012| , it will be shown that it is possible 
to considerably shorten the CPU time required to compute the 
diagram involved in the RegPT prescription with the use of a "fast" 
algorithm. The latter is based on the use of precomputed kernel 
functions making possible the computation of diagrams in ~ 0.01 
sees, per mode. We note that such a procedure can also be applied 
in the context of MPTbreeze. 



above a few h Mpc - . The use of these two-loop corrections 
improves upon the two-point propagator predicted value at 
least for z > 1 so that RegPT predictions can be potentially 
more precise. But it also makes the predictions less robust 
with regards to the overall spectrum as it makes the results 
more sensitive to nonlinear scales where baryon physics and 
shell crossings are likely to significantly affect the growth of 
structure. 



6 MULTI-POINT EXPANSION VS. N-BODY 
FOR VARYING COSMOLOGICAL MODELS 

In this section we extend the testing of the multi-point ex- 
pansion beyond the fiducial cosmological model used so far 
to measurements in our dedicated set of N-body simulations 
of various cosmologies described in Table [l] This kind of 
study is important given the number of assumptions leading 
to, say, Eqs. ( 29|30|31 l. It also helps setting in robust terms 
the validity and usefulness of our approach. Our results are 
summarized in Figs. [5] |6j [7] and [8] 

A clear picture of the different scenarios investigated is 
given in the top-left panel of Fig. [7] where we show the ratio 
of linear spectra of the simulated cosmologies to a reference 
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Figure 8. Power spectrum performance for different cosmologies at z = 1. Symbols correspond to simulation measurements, dashed red 
lines to linear theory, black solid to halof it and blue solid to the multi-point expansion, which remains accurate at < 2% on at least all 
scales showing Baryon Acoustic Oscillations. 
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Figure 9. Performance of the model presented in this paper for Las Damas cosmology (detailed in Table [TJl as a function of redshift. 
Data points with error bars are the measurements in three runs of Oriana. Solid blue are the predictions by the multi-point expansion 
(this work), while red dashed and solid black corresponds to linear theory and halof it respectively. 



smooth poweipj The amplitude of the longest wavelength 
modes (i.e. P x/ ' i ) varies by up to 70%, and a similar spread 
is found for the effective tilt at k ~ 0.01 /iMpc -1 . 

Let us first concentrate in the prediction for the two- 



simple model given in Eq. ( 16 1 performs at the sub-percent 



point propagator given in Eq. ( 16 1 for the different cosmolog- 



ical models and its comparison to simulation measurements. 
This is depicted in Fig.[5]for propagators measured at z = 
and Fig. [6] for those at z = 1. Remarkably in all cases the 



16 The smooth baseline used throughout the paper is a BBKS 
transfer function pardeen et al.||l986| l of shape F = 0.135 with 
tilt n s = 0.99 and arbitrary normalization. 



level for all the scales of interest, at least up to z = 1. 

We now turn to the analysis of power spectrum pre- 
diction using multi-point expansion for all the cosmological 
models detailed in Table [l] Figures [7] and [8] show the mea- 
sured power spectrum at z = and z = 1 respectively. Solid 
blue line is the two- loop model from Eqs. ( |29|31[ ), solid black 
corresponds to halof it and dashed red to linear theory. 

At both z = and z = 1 the multi-point expansion 
perform as for the FID case, that is, it matches the mea- 
surements up to k ~ a^ 1 (z) (for the given cosmology) 
at the < 2% level. This means, broadly speaking, up to 
k = 0.15/iMpc" 1 at z = and k = 0.25/iMpc" 1 at z = 1. 
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In turn halofit works at the < 6% level. In some cases, 
however, this departure appears at rather low k (e.g. for 
low-f2 m ). 

In Fig. [9] we concentrate in the performance of our im- 
plementation of the multi-point expansion as a function of 
redshift, from low z to high z using LasDamas measure- 
ments as a benchmark. Halofit seems to perform best 
around z ~ 1 but for smaller and higher redshifts departs 
more substantially from the measurements. For example, at 
z = Halofit is suppressed compared to simulations by 
5% at k = 0.15/iMpc -1 and 8% at k = 0.3/iMpc -1 . Us- 
ing Halofit as a benchmark for PT calculations is therefore 



not accurate enough for present-day work (Crocce & Scoc- 
cimarro|2008 1. It may be for this reason that Carrasco et al 
(20121 find a strong effect due to shell crossing at weakly 



nonlinear scales, after checking their effective stress tensor 
reproduces Halofit for the LasDamas cosmology. In fact, 



in earlier work Pueblas & Scoccimarro (2009) characterized 



the impact of a non-zero stress tensor and concluded that 
(extrapolating to the LasDamas cosmology at z = 0), shell 
crossing effects are about 1% (2%) at k = 0.2 (0.3) /iMpc" 1 



(see also |Pietroni et aL (2012 1). This justifies ignoring these 
effects in the calculations we present here at the scales of 



validity of MPTbreeze (note that as shown in Pueblas & 
Scoccimarro ( 2009 1 the effects of shell-crossing are strongly 



redshift-dependent, thus at z > they are less of a concern). 

It is worth stressing here that in all cases the perfor- 
mance of our code is maintained at a few seconds, as in the 
FID case described in Sec. 15.11 

Error bars shown in all power spectrum figures corre- 
spond to the expected statistical error (variance of the band 
power spectra) for the given simulation box-size. We have 
found that using only four realizations to estimate this error 
can be very unreliable. One needs at least ten or more runs 
to estimate this robustly (we have perform this testing using 
our 50 FID runs). Instead we found that the FKP expression 
( Feldman et al.|1994 l below matches the resulting ensemble 
error in our large FID ensemble very well on the scales we 
are interested (k < 0.4 h Mpc -1 ) 



crp 
P 



{Auk 2 5k) /k)' 



(34) 



where Sk is the particular binning used in the P estimation 
and kf = 27r/Lbox- Hence we use Eq. ( |34[ | to depict error 
estimates in Figsj7]j8]and[9] 

In Figs. |4l7|8l |9{ we have decided to show, besides 
MPTbreeze results, those from tools with comparable func- 
tionality such as halofit. We have nonetheless tested that 
RPT, as detailed in |Crocce fc S coccimarro (2008), is also ac- 
curate at the percent level but in an broader range of scales. 
However this is at the expense of longer evaluation times 
(hours as opposed to seconds). 

In addition we note that a fitting formulae for the non- 
linear power spectrum that updates halofit is available, the 
Coyote Universe emulator ( Lawrence et al.|2010 1 . However, 
this emulator can only be used within a specific region of pa- 
rameter space which unfortunately exclude all the cosmolo- 
gies investigated in our paper, even the one of LasDamas. 
This is mostly because the Coyote emulator does not treat 
h as an independent variable, but rather as one fixed by 
the distance to the large scattering surface as measured by 



CMB. Hence in Figs. [7] [8] and [9] we can only show the esti- 
mates from halofit. 



7 CONCLUSIONS 

With MPTbreeze we have implemented a renormalized per- 
turbative approach for the nonlinear power spectrum, the so 
called multi-point propagators expansion. We put particu- 
lar emphasis on the description of the BAO range of scales 
from low to high redshift. The main advantage over other 
techniques already present in the literature is that the eval- 
uation time is of the order of 5-10 seconds (in a single CPU), 
as discussed in Sec. 15.11 

Our implementation is based on a phenomenological 
description of the multi-point propagators themselves. In 
particular how their scale-dependence interpolates between 
their perturbation theory forms at low-A; and their large- 
k behavior obtained from non-perturbative re-summations. 
For the two-point (nonlinear) propagator at late times there 
is an unambiguous way to do this interpolation through well- 
known one-loop results, this is discussed in Sec. |4.1| 

We compared the adopted prescription with propaga- 
tor measurements in N-body simulations for seven differ- 
ent cosmological models (listed in Table [I] and discussed in 
more detail below). Remarkably we find it always gives a 
sub-percent agreement for all scales of interest, as shown in 
Figs, [l] [5] and [6] In addition, we developed simulations with 
independent initial positions and velocities. This allowed us 
to test for the first time the full matrix structure of the 
two-point propagator. We found that all individual compo- 
nents show an exponential suppression towards small scales, 
in turn very well reproduced by our simple prescription (see 
Fig.§. 

We then moved to the three-point propagator r*- 2 - 1 that 
is a function of triangle configurations. Hence its description 
is a priori much more complex. Recently |Bernardeau et al.| 
( 2012 1 put forward an interpolation scheme that respects the 



low-fe at arbitrary loop order and the high-fc limit. Evaluat- 
ing this would require one-loop calculations of r' 3 ', which 
we find slow to integrate numerically. Instead, we found that 
the decay of the three-point propagator when compared to 
measurements in N-body simulations is to a large extent the 
same as for the two-point case described above (in turn much 
faster to evaluate). In particular for configurations relevant 
for P(k) calculations, see Fig. [3] 

Provided with these results for the multi-point propa- 
gators we implemented the MP expansion in Eq. (|2| for the 
nonlinear power spectrum. In parallel we developed a set 
of large N-body simulations of different cosmological mod- 
els. Given the number of approximations in implementing 
a practical approach, testing against different cosmological 
models seems the only route to establish robust conclusions. 
The top-left panel of Fig. [7] gives an idea of the different 
slopes and amplitudes of the linear spectrum simulated in 
this paper. 

Comparison of P(k) measurements at various redshifts 
and cosmological models and our theory modeling is pre- 
sented in Figs. [7J [8] [9] Overall we can establish the accuracy 
of our approach at the 2% leve l on weakly nonlinear scales 
(roughly up to a^ 1 , given in Eq. 10 1. It always improves over 
standard PT and halofit. As mentioned above the evalu- 
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ation time scale remains around 5 to 10 sees, depending on 
redshift, integration accuracy, etc. 

Probably the main disadvantage of our implementation 
is that the range of validity does not extend much beyond 
BAO scales, particularly a low-z. This can be overcome 
by interpolating weakly nonlinear scales, as described by 
our approach, with high-fc asymptotic given by halo models 
(Vala geas &i Nishimichi||2011[) or P(k) resummation tech- 



niques ( |Anselmi fc Pietroni||2012 1. We leave this extra con- 
struction for future work. 

The code used to compute the multipoint expan- 
sion presented in this work is publicly available at 
http: //maia. ice . cat/crocce/mptbreeze/ 
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APPENDIX A: THE EDS APPROXIMATION 

The derivation of the multi-point expansion as presented 



Bernardeau et al. ( 2008 1 assumes that the cosmologi- 
cal model is such that f2 m (r)//+(r) 2 = 1, where f+(r) — 
d In D+ (r)/d In a and D+(t) is the growing mode solution 
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for the density contrast of the linearized equations of mo- 
tion: 5(k,r) = D+(r)5 (k). 

In this approximation the equations of motion simplify 
considerably reducing basically to the ones in an Einstein de 
Sitter background (fl m = 1) with a factorized linear growth 
factor, i.e. replacing the one in EdS, a(r), for the corre- 
sponding D+(t) (see Scoccimarro et al. (19981). Then at 



each perturbative order the solution becomes separable in r 
and k with the corresponding most growing term satisfying 
D n = (-D+) n and the PT kernels reducing to the ones in an 
EdS Universe. 

This approximation is known to be very accurate be- 
cause for most of the cosmic evolution fl m ~ 1. It is usu- 
ally followed in standard perturbation theory as well in re- 
summed approaches such as those in|Crocce fc Scoccimarro| 
( 2006b}, |Valageas| (|2007a[) |Matarrese fc Pietroni| ( |2007[ ), 
Taruya fc Hiramatsu| (|2008[> and |Matsubara| (|2008^ 



Pietroni ( 2008 1 investigated its limitations by numeri- 



cally integrating a system of coupled differential equations 



involving the power spectrum and bispectrum (see also Hi 



ramatsu & Taruya (2009 1). Here we follow an alternative 
approach implementing numerical simulations to directly 
address the validity of this approximation, in particular at 
BAO scales where we expect our analytical model to yield 
percent-level predictions. A different testing using numerical 



simulations is presented in McDonald et al. (2006) 



From the theoretical point of view the approach de- 
scribed above amounts to say that all the (nonlinear) cosmo- 
logical evolution is set by the linear growth factor D+ . Hence 
two cosmological models with the same linear spectrum 
(normalized at say z = 0) will show the same nonlinear P(k) 
at times n and r 2 such that £> +imo dell(*7"i) = D+,model 2(7-2). 

Hence we did the following numerical experiment: pro- 
vided with our fiducial LCDM run for fi m = 0.27 started at 
^lcdm _ -^y e implemented a complementary CDM sim- 
ulation (fl m — 1) started at a time that matched the growth 
from the initial conditions in the LCDM case. This yielded 
z com = 37 We choge outputs for the LCDM at z = 0, 0.3 

and 1 with corresponding growth from the initial redshift 
D+ = 38, 32.78 and 23.66. Thus we output the Q m = 1 CDM 
run at times that matched these growth factors, z = 0, 0.159 
and 0.605. We then compared the measured power spectrum 
and two-point propagator in the LCDM run with the cor- 
responding "matched growth" one in the CDM case. If the 
approximation were perfect these ratios would be 1. 

Results are shown in Fig. |Al| for the redshifts mentioned 
above. Top panel corresponds to the power spectrum com- 
parisons and bottom to the nonlinear two-point propagator. 

As expected on the largest scales (where linear the- 
ory applies) the ratio is indeed unity. At high z (where 
Q. m is closer to 1) the ratio is still unity within 0.2% for 
k < 0.4/iMpc -1 In turn, at low z the approximation does 
break down towards small scales. Nonetheless BAO scales 
are mostly unaffected. For instance, at k < 0.2 hT 1 Mpc and 
z = 0.3 (z = 0) the CDM power is < 0.5% (< 1%) smaller 



than the LCDM one. Pietroni ( 2008 1 finds similar devia- 



tions but a factor of two smaller (0.5% error at z = and 
k ~ 0.2/iMpc^ 1 ). 

In any case the approximation is found accurate at the 
sub-percent level on large BAO scales at low-2 with extended 
validity at higher z and smaller scales. 
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Figure Al. An approximate solution of the PT equations of 
motion (at each order) for an arbitrary ACDM background is 
obtained by solving the CDM case with a time dependence 
given by the linear growth of the ACDM model, i.e. changing 
a(r) — > D+(r). We test this approximation here by running 
"growth matched" simulations of similar ACDM and Q m = 1 
CDM models (except by their value of f2 m ). The top panel shows 
the relative difference in the measured power spectra and the bot- 
tom panel the ratio of nonlinear propagators. The approximation 
is weakest at low-2 but it is always never worse than 1% at BAO 
scales. 
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